

library(MASS)
library(tseries)

# pour les indices europeens
data(EuStockMarkets)
cac <- EuStockMarkets[,"CAC"]

# on affiche l'acf et la pacf de quelques series...
# on remarque que la periodicite s'observe sur l'acf

# pour afficher les diagnostics d'un fit : tsdiag

# un bruit blanc
bb = rnorm(n, 0, 1)
layout(matrix(1:3, 3, 1))
plot.ts(bb)
acf(bb)
pacf(bb)

# On choisit une periode de simulation tres longue, 
# comme ca on a une bonne estimation de l'acf et de la pacf
# On retrouve bien les resultats obtenus en cours...
n = 2000


# un modele MA


myplot = function(X) {
 layout(matrix(1:3, 3, 1))
 plot.ts(X)
 acf(X)
 pacf(X)	
}

ma.simu = arima.sim(n, model=list(ma = c(3.5, 4.5, 1.5)))
myplot(ma.simu)




# un modele AR
ar.simu = arima.sim(n, model=list(ar = c(0.9, -0.3) ))
layout(matrix(1:3, 3, 1))
plot.ts(ar.simu)
acf(ar.simu)
pacf(ar.simu)

# un autre modele AR
ar.simu = arima.sim(200, model=list(ar = c(0.9, -0.6) ))
layout(matrix(1:3, 3, 1))
plot.ts(ar.simu)
acf(ar.simu)
pacf(ar.simu)


# un modele ARMA

arma.simu = arima.sim(n, model=list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
          sd = sqrt(0.1796))
layout(matrix(1:3, 3, 1))
plot.ts(arma.simu)
acf(arma.simu)
pacf(arma.simu)




# un modele ARIMA 
arima.simu = arima.sim(n, model=list(ar=0.7))
myplot(arima.simu)
#arima.simu = arima.sim(n, model=list(order = c(2, 1, 2), ar=c(0.2))  )
layout(matrix(1:3, 3, 1))
plot.ts(arima.simu)
acf(arima.simu)
pacf(arima.simu)
# on constate une decroissance tres lente de l'ACF, car pas stationnaire

## INTERESSANT ICI


# proximite entre un ARIMA et AR

X = arima.sim(n, model=list(order = c(1, 1, 0), ar=0.8), sd=1)
myplot(X)


# the differenced series
dX = X - lag(X)
myplot(dX)

dX.fit = arima(dX, order=c(1, 0, 0))
# on retrouve bien les bonnes valeurs

# si on essaie de fiter un AR(2) a X
arima(X, order=c(2, 0, 0))
# on retrouve des coefficients qui correspondent a celui de l'ARIMA non stationnaire
# d'ailleurs on obtient un message d'erreur

# arima() ne fait pas la differentiation toute seule...



# un modele GARCH ?

layout(matrix(1:3, 3, 1))
plot(lh)
acf(lh, ci=0)
pacf(lh, ci=0)

layout(matrix(1:3, 3, 1))
plot(deaths)
acf(deaths, ci=0)
pacf(deaths)


# returns des indices europeens
data(EuStockMarkets)
cac = EuStockMarkets[,"CAC"]

cac <- log( EuStockMarkets[,"CAC"] / lag(EuStockMarkets[,"CAC"]))
dax <- EuStockMarkets[,"DAX"]
smi <- EuStockMarkets[,"SMI"]
ftse <- EuStockMarkets[,"FTSE"]

# SP500
SP500 # il s'agit des returns
idx = cac

# daily returns of SP500
layout(matrix(1:2, 2, 1))
sp <-  exp(cumsum(idx)/100)
ts.plot(exp(cumsum(idx)/100))
rsp <- idx
ts.plot(idx, main="Returns")

# ACF of SP500 and its returns
layout(matrix(1:3, 3, 1))
plot.ts(idx)
acf(idx, ci=0)
pacf(idx, ci=0)


layout(matrix(1:2, 2, 1))
acf(idx, ci=0, main="ACF of SP500 returns")
acf(idx^2, ci=0, main="ACF of SP500 returns^2")


# analyse spectrale, on peut trouver la frequence en regardant le periodogramme
# ou sa version lissee
layout(matrix(1:2, 2, 1))
myplot(deaths)
spectrum(deaths, span=c(3, 3))
spectrum(deaths)

## faire une stl, mais ca fait quoi au juste ???
plot(stl(deaths, 12))

# on fitte un AR sur lh
layout(matrix(1:3, 3, 1))
plot(lh)
acf(lh, ci=0)
pacf(lh, ci=0)

lh.fit = ar(lh, aic=F, order.max=1)
cpgram(lh.fit$resid, main="AR(1) fit for lh")

lh.fit = ar(lh, aic=T)
cpgram(lh.fit$resid, main="AR(3) fit for lh")
tsdiag(lh.fit)

lh.fit2 = arima(lh, order=c(3, 0, 0))
tsdiag(lh.fit2)


# prediction de la suite de la serie lh sur une autre periode de 12 mois
lh.fore = predict(lh.fit2, 12)
ts.plot(lh, lh.fore$pred, lh.fore$pred + 2*lh.fore$se, lh.fore$pred - 2*lh.fore$se )

